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Abstract 

The multigrid algorithm is a multilevel approach to 
accelerate the numerical solution of discretized dif- 
ferential equations in physical problems involving 
long-range interactions. Multiresolution analysis of 
wavelet theory provides an efficient representation 
of functions which exhibit localized bursts of short 
length-scale behavior. Applications such as comput- 
ing the electrostatic field in and around a molecule 
should benefit from both approaches. In this work, 
we demonstrate how a novel interpolating wavelet 
transform, which in itself is the synthesis of finite 
element analysis and wavelet theory, may be used 
as the mathematical bridge to connect the two ap- 
proaches. The result is a specialized multigrid algo- 
rithm which may be applied to problems expressed 
in wavelet bases. With this approach, interpolation 
and restriction operators and grids for the multigrid 
algorithm are predetermined by an interpolating mul- 
tiresolution analysis. We will present the new method 
and contrast its efficiency with standard wavelet and 
multigrid approaches. 



1 Introduction 

Many applications such as electronic structure calcu- 
lations require an efficient method to solve discretized 
linear differential equations which are generators of 
long range interactions, such as Poisson's equation. 
In modern electronic structure calculations, except 
those using planewaves, solving Poisson's equation 
takes up a significant amount of computational time. 
It is essential to produce an algorithm that has the 
capacity to solve efficiently problems involving such 
long range interactions. 

The multigrid algorithm is the state of the art in 



solving discretized linear differential equations de- 
scribing long range interactions. It uses a basic it- 
erative method, such as the weighted Jacobi method, 
over a sequence of scales to iterate faster toward the 
exact solution. In some cases, even multigrid does 
not produce sufficiently rapid convergence, particu- 
larly in calculations using wavelet bases. 

Multiresolution analysis, from wavelet theory, pro- 
vides the ability to carry out calculations on non- 
uniform grids, focusing resolution only in regions 
where it is needed. It is especially useful if the local 
region which requires high resolution moves through 
space during the calculation, as when the atoms move 
in a molecule in an electronic structure calculation. 
However, a basis set alone does not provide compu- 
tational efficiency. 

Often, a finite element basis set provides the ana- 
lytic framework for multigrid calculations. Combin- 
ing the advantages of finite element analysis and the 
multigrid algorithm with wavelet theory, could pro- 
vide a very powerful method to solve discretized lin- 
ear differential equations. This present research aims 
to unify these approaches for use in a wide range of 
applications, especially in electronic structure calcu- 
lations, where all of the aforementioned approaches 
prove to be very beneficial in various stages of the 
calculations individually. 

A basis set that combines wavelet theory with fi- 
nite element analysis exists and is known as the in- 
terpolet basis (Interpolets is based upon Deslauriers 
and Dubuc scaling functions [4], which play the role 
of both scaling functions and interpolets in the in- 
tcrpolet basis. See also [5]). The present research 
forms a natural bridge in the synthesis of the multi- 
grid algorithm with this basis, thereby adding the 
power of multiresolution analysis to the multigrid al- 
gorithm. In particular, we introduce this new Hop- 
grid algorithm for use in problems where an under- 



1 



lying wavelet or interpolet basis is necessitated by 
the problem at hand and application of multigrid al- 
gorithm could provide very useful in solving linear 
differential equations in the problem. 

In Section ^ of the paper, to establish a common 
notation, we give a brief overview of the traditional 
multigrid algorithm and the iterative methods that 
form the framework for this algorithm. Section ^ in- 
troduces the interpolet theory. Section ^ includes the 
description of non-orthogonal multiresolution analy- 
sis from wavelet theory in the interpolet basis. Sec- 
tion H presents the basic ideas behind expressing the 
multigrid algorithm in the interpolet basis, and de- 
scribes how such a union becomes more efficient by 
using multiresolution analysis. Section |^ introduces 
the new algorithm formed by a more impact synthesis 
of the multigrid algorithm and the interpolet basis, 
the Hopgrid Algorithm. The last section expresses 
the results obtained using this new algorithm. The 
appendix of this paper aims to familiarize the reader 
with the essentials of the multigrid algorithm with 
a thorough discussion of the basic iterative methods 
and the framework of this algorithm. 

2 Multigrid Algorithm 

This section briefly summarizes the multigrid algo- 
rithm [2] and the interpolation and restriction oper- 
ators used in this algorithm. A more detailed de- 
scription of the multigrid algorithm and these opera- 
tors, aimed for the physicist, may be found in Section 
^. Throughout this section, superscripts are used to 
refer to the scale on which a particular level of the 
problem is solved. 

2.1 Full V-cycle Multigrid Algorithm 

The multigrid algorithm is used to solve discretized 
linear differential equations of the form, 

Au ^ /, (1) 

where A is the linear operator, u is the solution vec- 
tor, and / is the source vector. Multigrid is based 
upon a basic iterative method, which is then applied 
over a sequence of scales in order to improve the con- 
vergence rate. One such iterative method used for 
this purpose is the weighted Jacobi method, which 
provides the following recursion, 

= W(„) + wD-^r^^n)- (2) 

Here, w is an appropriate weight (usually chosen to be 
w — 2/3), is the initial guess, D is the diagonal 



matrix composed of the diagonal elements of A, and 
r(„) is the residual at step n, defined as r(„) = f — 
Av(^n)- We refer to application of this recursion as 
relaxation. 

Multigrid exploits this basic iterative method in a 
multiscale fashion through the following procedure. 
First, a relaxations are performed on the scale on 
which the problem is defined. Then, the residual 
r^"^ is transfered up to a spatially coarser scale, us- 
ing a linear restriction operator, r^^'^ — R^^^^r''^^ On 
this coarser scale, one then solves the error equation 
^(1)5(1) = r^^\ relaxing again a times. (Note that, 
as discussed in the appendix, the exact solution to 
the error equation, when summed with the current 
iterative solution v, yields the exact solution to the 
problem u.) The procedure of passing up the residual 
to successively coarser scales continues until a prede- 
termined coarsest scale n is reached. At this coarsest 
scale n, one then transfers the solution vector t;'^"^ 
down to the next finer scale n — 1 using a linear inter- 
polation operator I^^^ . The error equation is then 

relaxed on this scale P times, using t;("^i) '^-^(n) 
as the initial guess. This procedure is followed down 
to the finest scale, relaxing (3 times at each scale i 
with the initial guess i;^*-* + I(j^^i-jV^^^^\ where finally 
one relaxes f3 times on the original linear equation to 
obtain an approximate solution to the overall prob- 
lem. 

The procedure described above is the V-cycle of 
the multigrid algorithm. V-cycles may be repeated in 
succession to produce a solution of any desired accu- 
racy. 

2.2 Theory for Interpolation and Re- 
striction Operators 

The interpolation and restriction operators should 
obey two conditions, known as the variational prop- 
erties. The first condition, the Galerkin condition, 
requires that 

^^(«+i)^(„).(n) ^ ^ 

(n) ^ ' 

It is a recipe for the linear operator to be used in 
relaxations on the next coarser scale. The second 
condition states that 

where c is a scalar constant. This is the recipe for the 
restriction operator, up to an overall scalar constant, 
given the interpolation operator. 
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3 Interpolet Theory 

In this section, we briefly review interpolet theory 
based on the discussion in Lippert, Arias and Edel- 
man [1]. 

In interpolet theory, functions varying slowly over 
integer length scales can be closely approximated as 
linear combinations of interpolating functions. 



(5) 



where the /„ are the expansion coefficients and the 
I{x — n) are functions with compact support, also 
known as interpolets. 

Interpolets are a basis set combining wavelet the- 
ory with finite element analysis. In place of the or- 
thonormality condition common in wavelet theory, 
interpolets have cardinality and interpolation from 
finite element analysis as their characteristic proper- 
ties. We will first discuss the properties which in- 
terpolets share with finite elements and then we will 
describe those properties which they share with tra- 
ditional wavelets. 

Cardinality means that the values of a function are 
zero at all integers except for zero, 



I{n) — 6no, for all integers n. 



(6) 



As a consequence of this condition, the function 
formed by the linear expansion in Eq. (||) will match 
exactly the value of the original continuous function 
at the integer grid points when the expansion coeffi- 
cients are taken to be the values of f{x) at the inte- 
gers, /„ = f{n), 

f{n)^^fm^{n-m)^^f„iSnm = fn- (7) 

m m 

Interpolation is the further condition that Eq. (^) 
reproduces any polynomial, up to order L for all x. 



^n''XL{x~k), /c = 0,1,...,L, 



(8) 



where subscript L denotes the order of interpolation. 
Given the values of a function at the integers, inter- 
polets then provide a compact estimate correct to L*^ 
order of the values of the function at the half integers 
through 



(9) 



The matrix representations for the interpolation op- 
erators defined by this relation for both first order in- 
terpolets {L = 1) and third order interpolets {L = 3) 



are given below. Note that the L = 1 case gives the 
familiar linear interpolation procedure, and that the 
i = 3 case is distinct from the interpolation given 
by traditional cubic i?— splines. Interpolet theory 
thereby provides guidance in selecting natural inter- 
polation operators for use in the multigrid algorithm. 
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In addition to properties from finite element analy- 
sis, interpolets satisfy the central condition of wavelet 
theory, the two-scale relation. This relation states 
that every coarse scale interpolet is expressible as a 
linear combination of finer scale interpolets. Interpo- 
lets therefore interpolate themselves exactly, 



Il{x) 



'c„Xi(2a; - n). 



(10) 



Here, the Tl{x) is the coarse scale interpolet, the 
XL{2x—n) are fine scale interpolets, and the c„ are the 
expansion coefficients. (Note that from the cardinal- 
ity property, we have that the c„ are just the values 
of the interpolet at the half-integers c„ = TL{n/2).) 
The significance of this condition is to ensure that a 
basis made from interpolets of varying scales always 
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Table I: Values of interpolets of orders ^ = 1, 3 at var- 
ious points X within their compact support. Values 
not shown are zero. 

o o o o 

Basis spanning the one dimensional space (Vo), whose resokiiion we wani io double. 

"'OOOOOOOO 

Basis with doubled spatial resolution spanning the space (VI), direct representation. 

.QoOoOoOo 

Basis with doubled spatial resolution spanning the space (Vo + Wo), MRA representation. 

■"0°o »0°o »0°o »0°o 

2 

Basis with resolution improved by a factor of 2 spanning a one dimensional space 

Figure 1: This figure illustrates pictorially how an 
MRA basis is created. The bigger circles represent 
the coarse scale interpolets, while the smaller circles 
are for the fine scale interpolets. 

provides a very uniform description of space. This 
point is described in detail the next section. 

The self-interpolation property may be applied re- 
cursively to express an interpolet in terms of interpo- 
lets on arbitrary finer scales, 

1l{x) = ... ^c„c,„...Cfc (11) 

n m k 

X ZL(2°2;-2"-in-2"-2^... -fc) 

= ^JL(fc/2°)JL(2°a;-fc). 

fe 

where a is a positive integer that determines the fine 
scale on which the original interpolet is expanded. 
This relation gives a procedure for determining the 
values of any interpolet recursively once given the 
values of the c„. Table ^ presents these expansion 
coefficients from which the interpolets may be con- 
structed through the two-scale relation. 

4 Non-orthogonal Multiresolu- 
tion Analysis 

In this section, we present a brief discussion of mul- 
tiresolution analysis (MRA) of interpolet and wavelet 



theories [1],[3], and the application of MRA to non- 
uniform grids such as the ones employed in electronic 
structure problems. 

Assume that there is a basis spanning a space Vq as 
in Figure |l|a. We refer to this space Vq as the coarse 
real space and this representation as the direct repre- 
sentation. In this space, any function /(x), varying 
slowly over the given grid spacing, can be approxi- 
mated with a linear combination of interpolets of the 
appropriate scale as explained in the previous section, 

f{x)=Y.^Mx-n). (12) 

n 

To double the resolution of this basis, we can double 
the number of grid points and decrease the scale of 
the basis functions by a factor of two, as Figure |^b 
illustrates. In this new space, Vi, which we refer to as 
the fine real space, any function f{x) varying slowly 
over the new grid spacing can be approximated with 
a linear combination of fine scale interpolets with half 
the compact support as 

/(x) = ^/„ZL(2x-n). (13) 

n 

General wavelet theory provides another tool to en- 
hance the spatial resolution of the basis in Figure |^a, 
locally. This tool is known as multiresolution analy- 
sis (MRA) . MRA states that to double the resolution 
of a given basis, we may add to the existing basis a 
new set of functions spanning a space Wq , such that 
Vq 4- 14^0 — Vi. For this condition to hold, we need 
to prove the two constraints that these spaces should 
satisfy, namely Vq-I-Wo ^ Vi and Vq-I-Wo CVi. In or- 
der for the first condition to hold, it is sufficient that 
the interpolets satisfy the two-scale relation, Eq. ( [l0| ) 
of Section |^. In order to prove the second condition, 
one must show that any interpolet Xl{2x — n) £ Vi 
can be written as a linear combination of interpolets 
in the space Vq -\- Wq. (This is proven for interpolets 
in [1].) The new basis spanning the space Vq -\- Wq is 
illustrated in Figure [Tjc. Henceforth, we will refer to 
such a mixed basis space as the MRA space and such 
a representation as the MRA representation. 

The MRA basis can be extended to include not 
only the next finer scale but other scales of increas- 
ing fineness, as Figure |I]d illustrates. Adding Wi^s 
{i = 0, 1, 2, TV) to Vq improves the resolution of 
the original space by a factor of 2^+^, 

Vn ^Vo + Wq + Wi + ... -\-Wi + ... + Wn-1- 

Multiresolution analysis is especially useful when 
the problem requires only local regions of high res- 
olution and only a small fraction of the expansion 
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coefficients for the MRA basis are significant. Un- 
der these circumstances, one needs to employ only a 
small subset of the full MRA basis. We shall refer 
to the points in space associated with the remaining 
basis functions as the non-uniform grid for the prob- 
lem. In electronic structure problems, for instance, 
only a small spherical region around the atomic core 
requires high resolution to describe the most rapid 
oscillations in the electronic wave functions. Such 
problems are best addressed using non-uniform grids 
where high densities of grid points are concentrated in 
concentric spheres with differing resolution surround- 
ing the atomic nuclei. The advantage of using such 
an MRA basis appears when the nuclei move. New 
grids do not have to be generated, and associated 
Pulley forces calculated, as the nuclei move. Rather, 
high resolution MRA basis functions may be simply 
turned on and then off as nuclei pass by. [1] provides 
a more detailed explanation of this application. 

5 Combining the Multigrid Al- 
gorithm with Interpolets and 
MRA 

We have so far described the multigrid algorithm and 
the interpolet basis. In this section, we discuss the 
rationale for combining the multigrid algorithm with 
interpolets. The result of Section |3.1| will be an al- 
gorithm in the direct representation. We will discuss 
multiresolution algorithms in Section BT2[ 



5.1 



Choice of Interpolation and Re- 
striction Operators 

Interpolets provide a natural prescription for an in- 
terpolation operator to be used in the multigrid algo- 
rithm. However, multigrid also requires an appropri- 
ate restriction operator so that both variational prop- 
erties are satisfied. In this section we will see that 
the appropriate restriction operator is just the trans- 
pose of the interpolet interpolation operator. This 
comes about because of the form of linear operators 
generated by applying the Galerkin technique to the 
interpolet basis. 

Let the following be the d-dimensional linear equa- 
tion to be solved. 



d<j>ix) = p{x) 



where O is a linear operator and (j){x) and p{x) are 
functions of the d-dimensional variable x. If we are 
to solve this equation in an interpolet basis, (/>(a;) and 



p{x) are then expanded as, 

4>{x) = y^^4>nlL{x - n), p{x) = '^PuTl^x - n). (15) 

n n 

Substituting these expansions into the linear equation 
and applying / dx II {x — m) to both sides yields 



(16) 



where 



A';^l = dxlL(x~ m) dlL{x - n) (17) 



is the interpolet form for the linear operator on scale 
n, and 

C/^"2 = [ dxlL{x- m)lL{x - n). (18) 



Using the two scale relation, we may relate this 
operator to the linear operator on the next coarser 
scale , 



dx - m)dlL{x/2 - n) 



CfcC; / dx 1l{x — 2m — k) OIl{x ~ 2n — I) 



2n+l 



Cl 



^(n+1) ^ 



(19) 



In the final line, we convert the relationship between 
the operators on the two scales into a matrix equa- 
tion. From Eq. ( [l9| ) , we may construct a set of opera- 
tors that automatically satisfy both variational prop- 
erties, the Galerkin condition and the full weighting 
condition. The appropriate set of operators are as 
follows. First, the linear operators ^^"^ representing 
O will be those generated from the interpolets using 
the standard Galerkin procedure. Second, the inter- 
polation operators II will be those generated from 
the interpolets as discussed in Section ^ Finally, this 
analysis shows that to complete the set, the restric- 
tion operators must be i? = ij. 

Implementing the interpolet basis in the multigrid 
algorithm is then using the linear operators A^"' to 
solve Eq. (|l]) and the error equation, and using the 
matrices and as the interpolation and the re- 
striction operators, respectively. 



(14) 5.2 Synthesis of MRA into the Inter- 
polet Multigrid Algorithm 



In this section, we discuss how we may generalize the 
algorithm that combines the interpolet basis with the 
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Original basis in the MRA representation. 
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Direct representation for the original basis. 

='00000000 

Direct representation of the restricted basis. 
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MRA representation of the restricted basis. 



Figure 2: Interpolation and restriction operations 
carried out in the MRA representation. The arrows 
indicate the direction to follow to carry out interpo- 
lation and restriction. The larger (smaller) circles 
represent the coarser (finer) scale interpolets. 



multigrid algorithm to multiresolution bases. We also 
discuss briefly how this multigrid algorithm in the 
MRA bases may be applied to non-uniform grids in 
an efficient manner. 

We will now discuss two approaches for carrying 
out interpolation and restriction for data in the MRA 
representation illustrated in Figure |^a. The first ap- 
proach is a direct approach and is equivalent to the 
method explained in the previous section. To carry 
out restriction in this approach, one first changes 
from the MRA representation to the direct represen- 
tation (the process carrying data from Figure [^a to 
Figure |^b. One then carries out the restriction op- 
eration in the usual way by applying the previously 
described linear operator (making the transfor- 
mation Figure |b->|c). Finally, one converts back 
from the direct representation to the MRA represen- 
tation (^— >|^). Interpolation within this approach 
also consists of three steps. First, we change from 
the MRA representation to the direct representation 
(^— >|^c), then we carry out interpolation with the 
usual operator II (§c^^)) and finally we convert 
back from the direct representation to the MRA rep- 
resentation (|^-^^). We refer to this approach of 
changing from the MRA representation to the direct 
representation in order to carry out the interpolation 
and the restriction operations as the operator method. 

The second approach is an improvement upon the 
operator method. This approach uses the properties 
of the MRA basis to carry out the operations of in- 
terpolation and restriction entirely within the MRA 
representation. We refer to this second method as the 
basis method. It uses the general physical meaning of 
restriction and interpolation to exploit the interpolet 



(o):operator method 
(x):basis method 



1 2 3 4 5 6 

# of V cycles 



(o):operator method 



(x):basis method 



WU (work units) 



Figure 3: The operator method compared to the basis 
method. The first figure represents the convergence 
rate vs the number of V-cycles and the second figure 
shows that the convergence rate of the basis method 
is computationally more efficient when the number 
of work units are considered. The above graphs are 
obtained by solving the Poisson's equation on a grid 
of 128 points, using first order interpolets. 
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basis. Restriction is the elimination of the fine in- 
formation from a given basis and mapping it onto a 
coarser one. In the interpolet basis, this fine infor- 
mation corresponding to the high frequency modes is 
carried by the finest interpolets. Therefore, restric- 
tion in this basis corresponds to the obviation of all 
the coefficients for the finest interpolets. Interpola- 
tion maps coarse information onto a finer grid. In 
the interpolet basis, the coarse information is repre- 
sented by the coarse interpolets. Thus, interpolation 
leaves the coarse interpolets unaffected, padding the 
coefficients for all finer interpolets with zeros. 

The operators in this new method satisfy the 
Galerkin condition since the interpolation operator in 
the MRA representation consists of a block identity 
matrix and a block zero matrix for the coarse and fine 
informations respectively. Figure ^ illustrates that 
the operator and the basis methods are equivalent 
with the same convergence rates per V-cycle. 

Note that the basis basis method has three distinct 
advantages over the operator method. First, the ba- 
sis method requires no matrix-vector multiplications 
to implement interpolation and restriction. Figure 
illustrates that the basis method is superior to the 
operator method on a per flop basis. Second, the ba- 
sis method provides the abilities to interpolate and 
to restrict over many scales in a single step. Finally, 
the basis method is easily generalized to problems 
on non-uniform grids. On non-uniform grids, restric- 
tion is merely obviation of finer scale coefficients and 
interpolation is padding finer scale coefficients with 
zeros. 

6 The Hopgrid Algorithm 

In this section we introduce the Hopgrid algorithm, 
which is the synthesis of the multigrid algorithm and 
the interpolet basis in the MRA representation. We 
present Hopgrid as an algorithm to solve discretized 
linear differential equations on interpolet bases in 
MRA representation. 

The efficiency of using the interpolet MRA basis 
in the multigrid algorithm appears when we look at 
how the error behaves at each step. The error vector 
at iteration step n is 

e(n)= A-^f -V(^n)- (20) 

Substituting the definition of the residual, r = f — 
Av and Eq. (||) into Eq. ( po|) returns the following 
expression, which states that error at each step is 
multiplied by a convergence operator 

e(n+i) = Ce(„) (21) 
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Figure 4: Eigenvalue spectrum of the convergence op- 
erators for the weighted Jacobi method in the tradi- 
tional multigrid basis and the interpolet MRA basis 
for a problem of 256 grid points. For the eigenvalue 
spectrum of the traditional weighted Jacobi operator, 
w=2/3, and for the weighted Jacobi operator in the 
interpolet MRA basis, w=l. 

where C = 1 — ■wD~^A. The eigenvalues of this op- 
erator determine how the error behaves at each iter- 
ative step. Figure ^ illustrates the eigenvalue spec- 
trum of the convergence operators for both the tra- 
ditional weighted Jacobi method and the weighted 
Jacobi method in the interpolet MRA basis. 

The convergence bands shown in Figure ^ are cho- 
sen for convenience to describe the effects of the 
eigenvalues on the error modes. The high convergence 
band is defined to be the region in which the choice of 
w for the weighted Jacobi recursion guarantees that 
about half of the eigenvalues of the convergence oper- 
ator of the traditional weighted Jacobi method lie in 
the vicinity of zero. The errors which are multiplied 
by these eigenvalues approach zero at a high rate. 
The medium convergence band and the slow conver- 
gence band are the regions where the convergence of 
the error to zero is slower. The distinction between 
these two bands in Figure ^ is chosen for purposes of 
illustration. 

First, we will consider the convergence of the tradi- 
tional weighted Jacobi method. The choice of weight 
w in the method affects the number of eigenvalues 
that lie in each of the bands; however, the operator 
D~^A always has eigenvalues approaching zero, and 
so the convergence operator always includes eigenval- 
ues near unity and modes which converge very slowly. 

On the other hand, most of the eigenvalues of the 
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Figure 5: Comparison between the weighted Jacobi 
method in the MRA representation and the Hopgrid 
algorithm for the solution of Poisson's equation on 
grids of several different sizes. 



convergence operator in the interpolet MRA basis are 
clustered in the high convergence band. The reason 
for this behavior of the eigenvalues is the fact that 
the flow of information on the interpolet MRA ba- 
sis is stronger than the flow of information in the 
basis which the traditional multigrid algorithm uses. 
The superiority in communication in the MRA ba- 
sis comes from the fact that the MRA basis includes 
functions with support spanning the entire spatial ex- 
tent of the problem. 

The information given in Figure ^ shows that af- 
ter one weighted Jacobi relaxation in the interpolet 
MRA basis, most of the error modes are eliminated. 
Only the modes that correspond to the eigenvalues 
at the edges of the spectrum survive after relaxation. 
Numerical experiments show that these surviving er- 
ror modes have low spatial frequencies and so may 
be eliminated very efflciently, using the basic ideas 
underlying the traditional multigrid algorithm. 

We now exploit these observations to produce a 
new algorithm. The basic outline of the algorithm 
is to first use weighted Jacobi in the MRA represen- 
tation to eliminate the error modes that correspond 
to the high frequency eigenvalues in the high conver- 
gence band, and then to restrict the residual to a very 
coarse spatial scale where the error equation is solved 
exactly, eliminating the remaining error modes. The 
fact that we hop over many scales in order to produce 
the error vector gives this new algorithm its name, 
the Hopgrid Algorithm. Below we shall refer to the 
scale on which we solve the error equation as the hop 



scale, h. Figure g| shows the comparison between the 
convergence rates in digits/WU (work unit) vs prob- 
lem size for the Hopgrid algorithm and the weighted 
Jacobi method in the MRA representation. As this 
figure illustrates, the hopping over many scales to 
calculate the error improves the convergence rate, es- 
pecially with increasing problem size. 

The detailed procedure for the Hopgrid algorithm 
is as follows. First, one applies weighted Jacobi in 
MRA representation once on the finest scale 
to find an approximate solution vector w^*^-*. Then, to 
reduce the computational overhead of the algorithm, 
instead of calculating all the entries of the residual 
7.(0) — J — Av^'^\ we calculate only the 2^ entries 
which correspond to the restricted residual on the 
hop scale h. Furthermore, because these entries de- 
pend only weakly on the high frequency components 
of v^^^ , we compute them only with the 2'*+^ entries 
of u'-^-' corresponding to the hop scale with two levels 
of refinement. After hopping over many scales, the 
number of grid points on h is very small compared 
to the number of points on the finest scale, so that, 
now, solving the error equation, ^('*)e('*) — r^'^^ ex- 
actly becomes a negligible part of the computation. 
Finally, we hop again over many scales and interpo- 
late the resulting e^''-' back to the finest scale, again 
just padding with zeros. The sum of the approximate 
solution before hopping v^^^ and the interpolated er- 
ror from scale n gives the improved solution vector 
V = -|-/^^|^je^"'. The hopping procedure may then 
be applied recursively to reach any desired level of 
accuracy. 

Numerical experiments show that the best conver- 
gence rates are obtained when the hop scale h is cho- 
sen so that /i < Y, where N is defined such that the 
number of grid points for the original problem is 2^. 

7 Results 

In this section we compare the efficiency of the Hop- 
grid algorithm with the traditional multigrid algo- 
rithm. We use two different measures of convergence 
rate for comparison. The first measure is the conven- 
tional measure in terms of digits per WU [2]. Here, we 
define one work unit to be the number of fiops it takes 
to multiply a vector by the matrix representing the 
linear operator in the corresponding basis. We will 
see that in terms of this measure the present algo- 
rithm is signficantly superior. The MRA representa- 
tion loses part of its advantage in problems requiring 
uniform resolution, as on a uniform grid, MRA op- 
erators are far denser, having fractal dimension 1.5, 
than their single scale counterparts. Even in this ex- 
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Figure 6: Source vector and the exact solution to 
Poisson's equation that is used in finding the results 
shown in this section. 

treme case, the Hopgrid algorithm overcomes the ap- 
pearant disadvantage of the density of the MRA op- 
erators. To demonstrate this we also compare the 
convergence rates of the two algorithms on a per flop 
basis. Note that we envision applying the Hopgrid al- 
gorithm to non-uniform problems expresses in MRA 
representations and the flop count for the uniform 
case represents a worse-case scenario. 

The results presented in this section all use third 
order interpolets. The differernatial equation we 
solve for our tests is the Poisson Equation in one di- 
mension on the unit interval with periodic boundary 
conditions, 

- -^cl^ix) = pix). (22) 

As Figure ^ illustrates, we use for testing purposes 
two (5-functions of weight one with opposite signs as 
the source vector p{x). Given this source vector, the 
exact (both analytic and numeric) solution to Pois- 
son's equation, (f>{x), is the piece-wise linear function 
as given in Figure 0. 

Following the procedure given in Section 5.1, we 
substitute the interpolet basis expansions for 0(x) 
and p{x) into Eq. (p^). In the direct representation, 
this procedure returns the matrix equation 

n n 

where A^nn = J dx ~ m)-£lL{x - n) and 

Umn = J dx Xl(x — m)TL (x ~ n). These are the ma- 
trix elements employed in the traditional multigrid 
calculations, numerical values computed by recursive 
application of the two scale relation as described in 
[1], are provided in Table ||. 



For Amn , values for \ m — n \ : 





1 


2 


3 


4 


5 


First order interpolet: 


-2 


1 











Third order interpolet: 


20 
9 


9 
8 





1 

72 









Table II: Non-zero matrix elements for the linear op- 
erator Amn for the first and the third order interpo- 
lets. The matrix elements for the fifth order inter- 
polets can be calculated in a straight forward manner 
but are not listed here. 




200 400 600 800 1000 1200 1400 1600 1800 2000 

Matrix size 



Figure 7: The standard matrix vector multiplication 
by the linear operator in the MRA basis vs the non- 
standard multiplication. 
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Figure 8: The convergence of the traditional multi- pjg^^^ g. Comparison between the convergence rates 
grid algorithm versus the convergence of the Hopgrid traditional multigrid algorithm and the Hop- 
algorithm for the solution of the Poisson's equation ^-^^ algorithm for increasing problem size, 
on a grid of 1024 points. 



10 



The Hopgrid calculations require matrix elements 
of the linear operator in the MRA represenation. 
These elements are 

AZn ^ I dx ^Il{2'^x - m)^lL{2'-x - n), (24) 
J ax ax 

where Si is the scale of basis function i. These values 
may be computed from the previous A^n by appro- 
priate change of integration variable and application 
of the two scale relation. 

Multiplying a vector by this MRA matrix may be 
performed using one of two different methods. The 
first is straight forward multiplication by the matrix 
j^mra^ Bccausc of the density of the MRA matrix 
mentioned above, the number of flops required to 
apply the operator within the method grows faster 
than linearly with the size of the problem. The 
second approach uses the two-scale relation to re- 
duce the number of operations so that the effort in 
the matrix multiplication scales only linearly with 
the size of the problem [1]. Figure |^ compares the 
number of M flops needed by the standard and the 
non-standard multiplications as a function of prob- 
lem size. We use this superior approach in all of our 
comparisons. 

Figure || shows how the convergence of the Hop- 
grid algorithm to the exact solution compares with 
the convergence of the traditional multigrid algorithm 
for a problem of 1024 grid points. In Figure ||a, the 
convergence to the exact solution is plotted against 
work units (WU). Using this comparison the Hop- 
grid proceesure is superior to the traditional multi- 
grid algorithm. When the convergence rates of these 
two algorithms are compared in terms of digits per 
Mflops, as in Figure ||b, we find that despite the in- 
creased density of the MRA matrix (by about a factor 
of ten), the two methods give comparable results on 
a uniform grid. 

Next we investigate the dependence of these con- 
vergence rates on problem size. In terms of digits 
per WU, Figure ||a shows that although the conver- 
gence rate of the traditional multigrid algorithm re- 
mains nearly constant, the convergence rate of the 
Hopgrid algorithm increases with increasing problem 
size. Figure ^ shows that even for the uniform case, 
the efficiency of the Hopgrid algorithm overcomes the 
increased complexity of the MRA matrix and results 
in nearly the same convergence rate as the traditional 
method in terms of a direct floating point operations 
comparison. It is in the case of non-uniform grids, 
where the application of the MRA matrix requires 
far fewer operations, where we expect the maximum 
from the new algorithm. 



8 Conclusion 

In this research, our aim was to develop a new algo- 
rithm that benefits from the best combination of a va- 
riety of approaches, in particular the multigrid algo- 
rithm, multiresolution analysis from wavelet theory, 
and finite element analysis. The application for the 
new algorithm which we have in mind is the numer- 
ical solution of problems involving long-range inter- 
actions where an underlying interpolet basis with its 
multiresolution properties is needed. We have shown 
that the interpolet theory, based upon multiresolu- 
tion analysis from wavelet theory and interpolation 
properties from finite element theory, provides a nat- 
ural and unique choice of interpolation and restriction 
operators and an underlying grid for the multigrid 
algorithm. Then we showed how the operations of 
restriction and interpolation may be carried out ef- 
ficiently, without the application of linear operators, 
by exploiting the properties of an MRA basis. Fi- 
nally, we introduced the Hopgrid algorithm, which 
draws upon these results to produce a multigrid-like 
method to solve discretized linear equations in prob- 
lems expressed in multiresolution wavelet bases. We 
have seen that the convergence rate of the new algo- 
rithm in terms of work units is superior to that of the 
traditional multigrid apprach and that, even in the 
worse-case scenario of a uniform problem, the new al- 
gorithm performs as well as the traditional approach 
on a direct floating point operation basis despite the 
increase complexity of MRA matrices. 

9 Appendix: Description of the 
Multigrid Algorithm 

In this section, we give a detailed discussion of the 
traditional multigrid algorithm. This introduction to 
multigrid is a review of the technique as given in A 
Multigrid Tutorial by William Briggs [2] . 

9.1 Basic Iterative Methods And the 
Coarse Grid Correction Scheme 

We first present the weighted Jacobi basic iterative 
method and the coarse grid correction scheme which 
are the core ideas in the multigrid algorithm. 

Assume that the discretized differential equation 
we wish to solve may be represented in matrix form 
as 

Au ^ /, (25) 

where A is the linear operator, u is the solution vec- 
tor, and / is the source vector. Given an approxima- 
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tion to the solution, which we will refer to as v, the 
error can be expressed as e — u~v. Another measure 
of error is the residual r, 



Eigenvalues of Convergence Operator 



slow convergence 



r = f — Av. 
With this definition, 

Ae = r 



(26) 



(27) 



is equivalent to Eq. (^5|). Henceforth, we shall refer 
to Eq. ( |2^ ) as the error equation. 

Decomposing the initial matrix A as A = D—L—U, 
where D, — i, —U are the diagonal, lower triangular 
and upper triangular elements of A respectively, Eq. 
(p5|) becomes 



medium convergence 




mode number 



D-HL + U)u + D-^f. 



(28) 



This equation may be solved iteratively with the fol- 
lowing recursion. 



Figure 10: Convergence operator for the weighted Ja- 
cobi method. 



(29) 



This recursive algorithm can be improved by choosing 
an appropriate weight w, and taking a weighted av- 
erage of the initial guess, «(„), and the solution from 
Eq. (H), 

W(„+i) = [(1 - w)I + wD-\L + [/)]«(„) + wD~ 

Solving Eq. (^6]) for / and substituting into Eq. (|30| ) 
gives a simplified form for this recursion, in terms of 
the residual alone, 



wD ^r 



(")' 



(31) 



where r(„) is the residual at step n. This simplified 
recursion is the core of the basic iterative method 
known as the weighted Jacobi method. The applica- 
tion of this method for solving a problem is referred 
to as relaxation. 

The weighted Jacobi method and similar basic it- 
erative methods rapidly eliminate high frequency er- 
rors while leaving the smooth components of the error 
mostly unaffected. This is called the smoothing prop- 
erty. The reason for this behavior can be explained 
when we look at the eigenvalues of the convergence 
operator C = 1 — wD~^A. The eigenvalues of this 
operator determine how fast the error converges to 
zero. Fig. |l^ shows that the eigenvalues of the con- 
vergence operator for the weighted jacobi method lie 
in various ranges. The bands are chosen for illustra- 
tion purposes and an explanation for these choices 
are given in the main body of the paper. 

The errors that are multiplied by the eigenvalues in 
the slow convergence band correspond to eigenvectors 



varying slowly in space. Because the convergence op- 
erator C does not eliminate efficiently most of these 
low frequency errors, this basic iterative method ap- 
proach the exact solution very slowly. A more effi- 
cient method would eliminate all the error modes at 
equal rates. The multigrid algorithm has this feature 
^ and produces a solution that approaches the analytic 

/. (30)gQiy^jQj^ with a higher convergence rate. 
1^-3. The transition from the weighted Jacobi method to 
the multigrid algorithm is made by improving upon 
the basic Jacobi method. One way to improve the effi- 
ciency of weighted Jacobi relaxations is to start with a 
better initial guess. Conceptually, one can do this by 
solving the problem on a spatially coarser scale and 
using this solution as the initial guess for the origi- 
nal problem. This procedure is referred to as coarse 
grid correction (CGC). On the coarse scale, there are 
fewer grid points, and on this grid the smooth er- 
ror modes from the fine scale now appear higher fre- 
quency modes. Replacing the original problem with a 
smaller one of mostly high frequency errors allows the 
application of the weighted Jacobi method to solve 
the problem more efficiently. This way, the conver- 
gence rate for the total error is accelerated since both 
the high and the low frequency error modes are elim- 
inated on the fine and the coarse scales. 

Figure |ri| illustrates the accelerated convergence 
rate of coarse grid correction scheme compared to the 
weighted Jacobi method. 



9.2 Full V-cycle Multigrid Algorithm 

The multigrid algorithm simply extends the idea of 
CGC over a sequence of scales. Throughout this sec- 
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(o):weighted Jacobi 
(•):coarse grid correction 



# of iterations 



(o);weigttted Jacobi 
{'):coarse grid correction 



0.005 0.01 0.015 



0.02 0.025 0.03 
# of IVlfiops 



0.035 0.04 0.045 0.05 



Figure 11: This graph shows the difference be- 
tween the efficiencies of a pure basic iterative method 
(weighted jacobi) and the coarse grid correction in 
finding the sohition for the Poisson's equation on a 
grid of 128 points. 



tion, superscripts are used to refer to the scale on 
which that particular level of the problem is solved. 

In order to find an approximation to the solution 
of the linear equation Au — f, we first apply the 
weighted Jacobi method as in Eq. ^ a times (do a 
relaxations) on the finest scale, to obtain an approx- 
imate solution w^^) and a residual r^"^ . The residual 
calculated on this scale is then transferred to the next 
coarser scale using a linear operator i?|oj, known as 

the restriction operator. In general, R^^^^ takes a 
vector from a fine scale i, and transfers it to the next 
coarse scale i + 1. (This operator is discussed later 
in the appendix.) For the restriction step, we assume 
that the relaxation on the finest scale has eliminated 
most of the high frequency error components, which 
lie in the null-space of restriction, so that little in- 
formation is lost in the restriction process. The re- 
stricted residual is now 



^(1) - 7?(i)^(o) 
Next we solve the residual equation 



(32) 



(33) 



where A^^^ is the linear operator on the coarse scale, 
performing a relaxations, using e*-^^ = as an initial 
guess. This relaxation produces an approximation to 
the error on the fine scale. We could add this er- 
ror to v^^^ to obtain an approximate solution for the 
fine scale problem, since u = v + e, as previously 
explained. However, the multigrid algorithm contin- 
ues to go down to coarser scales and uses CGC to 
find solutions for the error equations at each coarser 
scale to improve the convergence rate. Following this 
prescription, we compute the residual on scale 1 and 
transfer it to the next coarser scale, using the restric- 
tion procedure described. Repetition of such relax- 
ation and restriction steps continues until the pre- 
determined coarsest scale n is reached. During each 
step leading up to the coarsest scale, the generated 
approximate solution vectors, the v^^\ and the resid- 
uals, the r^^\ are stored to be used later in the algo- 
rithm. Although the notation used for the result of 
the relaxations at each scale is v, it is important to 
realize that except for the finest scale on which the 
original problem is defined, all the vectors w''-' are 
approximations to the error since the error equation 
is used during relaxations. 

At the coarsest scale, another linear operator, an 



(n-l) 



interpolation operator /^^^ 
imate solution vector w^") 



transfers the approx 
down to the next 



scale. In general, the interpolation operator /, 



finer 

(i+i)' 
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takes a vector from a coarse scale z + 1, and trans- 
fers it to the next fine scale i. (This operator is 
also discussed later.) The interpolated approximate 
solution vector, I^^^ is the approximation to 

the error for scale n — I. The sum of the approxi- 
mate solution and give the approx- 
imate solution on scale n — 1. However, the inter- 
polation process may introduce high frequency errors 
into this solution. We eliminate these errors by relax- 
ing P times on the residual equation using the initial 
guess ul" = + "'"'w'-"'' and the residual 
7.("-i) stored from the previous restriction steps as 
the source vector. This relaxation produces the final 
^("-1) which is then interpolated to the next finer 
scale. This procedure is followed down to the finest 
scale. On the finest scale, we relax /3 times on the 
original equation, Eq. (p5|), using the initial guess 
wl"-' = + lj-^}^v^^'' and / as the source vector. The 
solution obtained at this final step is the result of 
one iteration of the V-cycle of multigrid algorithm 
and gives an approximate solution to the problem. 

It is a common procedure to carry out only one 
relaxation per each scale. All the results presented 
here are done with a — (3 = \, and w — 2/3. 

9.3 Interpolation and Restriction Op- 
erators 

The interpolation and the restriction operators in the 
multigrid algorithm obey two constraints, known as 
the variational properties. The first of these con- 
straints is the Galerkin condition, which specifies the 
form of the linear operator on a coarser scale. As- 
sume that the error vector on scale n is in the range 
of interpolation, e*^"^ — /(^^]^-)e'"+^-', for some vector 
g("+i) on scale n -f 1. The residual equation at scale 
n, A(")e(") = r("), then reads 



/!(")/(") f,{ri+l) ^ (n) 
(n+l) 

Restricting both sides, we find 

(n) (n) ' 



or 



(34) 



(35) 



( ^6| ) with A^^^^^ as the linear operator and the re- 
stricted residual r'^'^~^'^\ we determine the error vector 
on scale n + l which when interpolated becomes the 
error vector we seek, e'-"-'. 

The second variational property specifies the re- 
striction operator upto a scalar constant, given the 
interpolation operator, — c{I^^^^^^)'^ , where c 

is a scalar constant. This condition gives the restric- 
tion operator the full weighting property. The full 
weighting restriction involves taking some weighted 
average of values with neighboring points, instead of 
simply down sampling, a process known as injection. 
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^("-Hi)g("-i-i) 



„(»+!) 



(36) 



where 



^(„+l) ^ ^(n+l)^(„) 

(ri) (n+l) 



(37) 



This is precisely the Galerkin condition for the lin- 
ear operator on the next coarse scale. By solving Eq. 
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